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ABSTRACT 

We present a probabilistic approach for inferring the parameters of the present day power-law stellar 
mass function (MF) of a resolved young star cluster. This technique (a) fully exploits the information 
content of a given dataset; (b) can account for observational uncertainties in a straightforward way; 
(c) assigns meaningful uncertainties to the inferred parameters; (d) avoids the pitfalls associated with 
binning data; and (e) can be applied to virtually any resolved young cluster, laying the groundwork 
for a systematic study of the high mass stellar MF (M > 1 M Q ). Using simulated clusters and 
Markov chain Monte Carlo sampling of the probability distribution functions, we show that estimates 
of the MF slope, a, are unbiased and that the uncertainty, Aa, depends primarily on the number 
of observed stars and on the range of stellar masses they span, assuming that the uncertainties on 
individual masses and the completeness are both well-characterized. Using idealized mock data, we 
compute the theoretical precision, i.e., lower limits, on a, and provide an analytic approximation for 
Aa as a function of the observed number of stars and mass range. Comparison with literature studies 
shows that ~ 3/4 of quoted uncertainties are smaller than the theoretical lower limit. By correcting 
these uncertainties to the theoretical lower limits, we find the literature studies yield (a) = 2.46, with 
a l-a dispersion of 0.35 dex. We verify that it is impossible for a power-law MF to obtain meaningful 
constraints on the upper mass limit of the IMF, beyond the lower bound of the most massive star 
actually observed. We show that avoiding substantial biases in the MF slope requires: (1) including 
the MF as a prior when deriving individual stellar mass estimates; (2) modeling the uncertainties in the 
individual stellar masses; and (3) fully characterizing and then explicitly modeling the completeness 
for stars of a given mass. The precision on MF slope recovery in this paper are lower limits, as we 
do not explicitly consider all possible sources of uncertainty, including dynamical effects (e.g., mass 
segregation), unresolved binaries, and non-coeval populations. We briefly discuss how each of these 
effects can be incorporated into extensions of the present framework. Finally, we emphasize that the 
technique and lessons learned are applicable to more general problems involving power-law fitting. 
Keywords: stars: luminosity function, mass function - galaxies: star clusters: general methods: 
statistical 



1. INTRODUCTION 

The high mass end of the stellar initial mass func- 
tion (IMF) underpins much of extragalactic astrophysics. 
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Stars more massive than a few solar masses are largely 
responsible for most chemical enrichment, dominate the 
spectral energy distribution blueward of ~ 1/xm for all 
star-forming galaxies, and are presumed to dominate 
stellar feedback processes on galactic scales. Conse- 
quently, the exact numbers and mass distributions of 
high mass stars are central to the interpretation of in- 
tegrated light from distant galaxies, chemical evolution 
models, the frequency of core-collapse supernovae, the 
evolution of star formation rates over cosmic time, and 
the efficiency of star formation on galactic sca l es (e.g. 



Leitherer et al.|[l999l [Bruzual fc Ch ariot 200 31 |Tinsley| 
19801 |Smartt||2009 i IMadau et al.||1996| |Schmidt||1959 



Kennicutt||1989[ |Leroy et al.||2008 ) . Further, the form of 



the IMF and its potential sensitivity to local environment 
carries significant implications for how stars form from 
dense molecular cores and how the individual stars affect 



the properties of stellar clusters (e.g., Elmegreen||1999 


2004; Wcidner & Kroupa 


2005 


Zinnecker &: Yorke||2007 


McKee & Ostrikcr 


2007; 


Wcicl 


ner et al. 2010 Portegies 


Zwart et al.||2010 


Bastian et ai.||2010p. 



of ways (e.g., Salpeter 1955 Miller fc Scalo 1979 Scalo 
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1986 Chabrier 2003|, in parti cular as a set of power- 
laws. Following Kroupa (2001), we adopt the following 
parameterization of the IMF: 



$(M) 



,dN, 
"dM 1 



d M~ 



MnA < M < Mh 



(1) 



where Cj are chosen to ensure continuity, $(M) is nor- 
malized such that 



J $(M) dM = 1 M , 



(2) 



and a = cui, « = 1, 3 is slope of each power-law component 
within specified mass interval such that 



ai = 0.3, 0.01 < M/Mq < 0.08 
a 2 = 1.3, 0.08 < M/Mq < 0.5 
a 3 = 2.3, 0.5 < M/Mq < M max , 

with Mmax being the maximal mass of a star that could 
have been formed. The value of M max can either be set 
by the maximal mass at which the stellar lifetime exceeds 
the (finite) age of the cluste r, or by the star- or cluster- 
formation process itself (e.g., Elmegreen |2^1|2004ljy^ri 
dner fe Kroupaj 12005} IGoodwin £ Kouwennoven||2U09| 



Myers et al.||20ll| |Krumholz et al.||2010[ |20lT^ The 
high mass portion of the IMF (i.e., a^) sets essentially 
all observables (e.g., luminosity, color, chemical enrich- 
ment, etc.) in extragalactic contexts, making analysis of 
a single power-law a reasonable approximation in most 
environments outside the Galaxy. We will therefore used 
shorthand notation a = 013 throughout this paper. 

Despite its widespread importance, the IMF for stars 
above ~ 1 M remains insufficiently constrained by ob- 
servations. As shown in Figure [l] and listed in Table 
[TJ measurements of the IMF slope from studies of re- 
solved star clusters exhibit a scatter of ~ 0.5 dex. Ad- 
ditionally, various compilation 'a-plots' from the liter- 
ature derive mean IMF slopes with vari ations of ~ 0.4 



Marchi et al. 


2010 


Bastian et al. 


2010) 



ments and variations in the average values of different 
compilations introduce significant syste matics into inter- 
preting the higher redshift universe (e.g., Cool et al. 2008 
Conroy et al.|20 09; Narayanan & Dave 2012]), and further 



inhibit differentiation among st ar formation models that 
result in a universa l IMF (e.g., |McKee fc Ostriker||2007 



vironmental conditions (e.g., Weidner & Kroupa 2005), 



each of which presen t a drastically di fferent v iew of how 



galaxies evolve (e.g., 


Pflamm-Altenburg & Kroupa 2009 


Haas & Anders 


2010 


). 



precision and accuracy of IMF slope measurements 
should be taken at face value, which either implies a 
formal rejection of a universal high mass IMF slope 
or whether it should be interpreted as a reflection of 
stochasticity in the measurements. 

Currently, the majority of high mass IMF studies di- 
rectly count stars or measure luminosity functions in sin- 
gle or small sets of clusters, with each employing different 




Figure 1. Selected literature measurements for the high mass 
(M > 1 Mq) IMF slope from resolved star counts (see Table IT}. 
Points have been placed on the x-axis at the mean mass of tne 
range covered by a given study. The thin grey lines indicate the 
span of the full mass range. Error bars in the y-direction reflect the 
quoted 1-cr uncertainties. The heterogeneity in the MF recovery 
techniques makes it difficult to discern whether the ~ 0.5 dex of 
scatter of the data is due to physical or observational challenges 
and further convolutes the significance of the uncertainties. 



observational and IMF recovery techniques. Many stud- 
ies selectively address critical issues including corrections 
for binary stars, spatially dependent observational com- 
pleteness, and uncertainties in individual stellar masses. 
Consequently, whether the broad range of measurements 
in Figure [T] are due to intrinsic physical effects or arc 
the result of systematics re mains an open and important 
question (see discussions in|Scalo||1986[ |Elmegreen||1999[ 

lmegreen|20091 IBastian et al.| 



Kroupa|2002l |Scalo|2005[ [El megreen|2009| |Bastian et 



2010) 



Myers et al.||20lT| ) and those that are sensitive to en- tir 



In practice, we can only measure the present day mass 
function (MF) of evolved clusters, in which the observ- 
able stellar masses are limited at the high mass end by 
the turnoff mass corresponding to the age of the clus- 
ter, although the object of ultimate importance is the 
entire IMF. For a simple stellar population, i.e., a star 
cluster, th e measured MF slope is i dentical to the IMF 
slope (e.g., Elmegreen fc Scalo|2006 1 modulo certain dy- 
namical (e.g., mass segregation) and observational (e.g., 
completeness) effects. Consequently, our understanding 
of the IMF critically depends on our ability to accurately 
and precisely constrain the MF. 

Systematically counting the relative numbers of in- 
dividual stars (with M > 1 Mq) in a large collection 
of young clusters has long been advocated as the op- 
timal approach to constraining the high mass MF (e.g., 
Scalo|1986||Kroupal2002"l|Scalo|2005l IBastian et al.|2010l 
Kroupa et al.||201ip. As improved datasets of resolved 



stars have emerged (e.g., |Zaritsky et al.||1997||Massey et 



al1|2006l IHoItzman et af.||^00ti| |Dalcant,on et ai-IIMf 



2012| |Bianchi et al.||2012^ , it is crucial that the analysis 
tools for constraining the high mass MF are in place to 
fully exploit the information content of the data and to 
properly account for all known uncertainties. 

Broadly speaking, the practical procedure of obtaining 
a MF constraint for a resolved young cluster consists of 
three steps. First, one obtains individual mass estimates 
for the stars, either based on photometry in conjunction 
with isochrones, or from spectroscopy. Second, one de- 
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termines the completeness of the sample for which mass 
estimates were obtained. Finally, one fits a MF model 
(e.g., Equation[T]) to the data to obtain confidence limits 
on the slope parameter, a. In its usual implementation, 
this approach uses % 2 minimization of a power-law model 
fit to a binned or cumulative representation of the mass 
function, N(M). The fit is taken over the masses deter- 
mined to be sufficiently complete and the error bars on 
each mass bin reflect Poisson statistics of the number of 
stars in the bin. 

However, there are severe limitations in this tradi- 
tional approach. Correlations between the number of 
stars per bin and the adopted bin size can introduce sig- 
nificant biases into the resulting MF slope and associ 



ated uncertainties (e .g., Mafz Apellaniz & Ubeda 



2005 



Cara & Lister 2008). Further, such an approach does 
not provide an objective framework for folding in obser- 
vational effects such as completeness and mass uncer- 
tainties, which are critical to accurate MF constraints. 
More robust approaches such as Monte-Carlo simula- 
tions, color-magnitude diagram fitting, non-parametric 
kerne l estimators, and maximum likelihood techniques 
(e.g. |Elson et al I|l989l |Gennaro et al.||2011| polphin 
2002 |Vio et al.|1994[|larrab|1982| ) alleviate some ot these 
issues. However, their reliability has not been shown in 
large, diverse sets of clusters, the interpretation of the er- 
ror bars is often ambiguous, and there remain statistical 
short-comings in even the most advanced conventional 
approaches (e.g., n oise amplification in maximum likeli- 



hood estimations; Gould et al.||1997j ). In short, tradi 



tional approaches to MF reconstruction are inadequate 
for a systematic investigation of the high mass stellar 
IMF. 

The main goal of this paper is to lay out a forward 
modeling technique that remedies at least some of the 
shortcomings of the traditional approach. It should re- 
sult in constraints on the high mass MF that come close 
to exploiting the full information content of the data and 
that explicitly account for known observational uncer- 
tainties. In particular, by laying out a forward model- 
ing formalism, we aim to overcome the following limita- 
tions of many previous analyses by: (1) accounting for 
the mass estimate errors of individual stars that may 
arise from the conversion of observed fluxes into stellar 
masses; (2) avoiding binning the measurement in mass, 
as the most massive bins are unavoidably sparsely popu- 
lated; (3) dealing explicitly with completeness functions 
that may not be a simple step function of stellar mass, 
and instead be either gradual or varying as a function of 
position; (4) implementing a rigorous way to derive the 
joint posterior distribution function of the MF parame- 
ters (e.g., a and M max ) in light of evidence for the data 
and its uncertainties and possibly other prior informa- 
tion. 

In this paper, we present the probabilistic framework 
for inferring the parameters of the high mass MF of a 
young stellar cluster. We carefully walk though the nec- 
essary mathematics and interlace illustrative examples 
using simulated clusters. For clarity in the derivations 
and examples, we have made several simplifying assump- 
tions. In particular, dynamical effects such as mass seg- 
regation, ejection, and relaxation are not included in the 
present paper, although their role in interpreting the ob- 
served MF as the IMF is indisputably important. 



The methodology presented in this paper is motivated 
in part by the Panchr omatic Hubble Andro meda Trea- 
sury program (PHAT; [Dalcanton et al.|2012[ ). This Hub- 
ble Space Telescope multi-cycle treasury program is map- 
ping ~ 1/4 of M31's star forming disk, providing near-UV 
through near-IR imaging of ~ 10 s resolved stars. Based 
on analysis of the Year 1 data, we anticipate the survey 
will resolve > 1000 young clusters (i.e., ages < 10 Myr; 
or a main sequ ence upper mass limit of > 5 M ; John- 
son et aL]|2012 ) with sensitivity down to stars with M ~ 
1-3 M , which will enable a large scale study of the high 
mass stellar MF. Although the terminology 'high mass' 
is often reserved for stars with M > 10-15 M Q , common 
IMF parameterizations (cf. Equation [lj suggest that all 
stars above ~ 1 M are drawn from the same underlying 
IMF, indicating that conventionally defined intermedi- 
ate and massive stars can be used to learn about their 
common IMF. Throughout the paper, we use these num- 
bers as guides for our simulations, but emphasize that 
the methodology developed in this paper can readily be 
generalized to measure the MF of any resolved cluster 
and for any parameterized MF model. 

The remainder of this paper is organized as follows. 
In Ej2j we lay out the general probabilistic framework for 
inferring the MF of a young stellar cluster. We then ap- 
ply this technique to highly idealized mock data (e.g., 
no mass uncertainties, perfectly known com plete ness) 
and present illustrative results in £j3] In §3.2| and §3.3[ 
we derive the theoretical lower limits for recovered MF 
slope precision and compare the results to select litera- 
ture studies. We then consider the more complex case 
of masses with finite mass uncertainties in <|4j and ex- 
plore the case where th e ste llar mass distributions are 
log-normal functions in §4.1| We then consider the case 
of linear completeness functions in f|5] Finally, we dis- 
cuss several caveats to the MF models adopted in this 
paper (e.g., the influence of cluster dynamics, non-coeval 
populations, unresolved binaries) in [J6j 

2. PROBABILISTIC FRAMEWORK FOR MODELING THE 
STELLAR MASS FUNCTION OF A CLUSTER 

We begin by illustrating the scope of the problem at 
hand. Consider the IMF of a hypothetical cluster as 
shown in Figure [2] We presume that the high mass IMF 
regime of this cluster consists of some number of stars, 
whose mass distribution follows a power-law with a slope 
of a, and whose masses are bounded by M mm and M max , 
which is the mass of the most massive star that could 
have formed in the cluster, and may not be the same for 
different clusters. This 'zero-age' state of the cluster is 
the desired dataset to make direct inferences about its 
stellar IMF. 

However, there are several considerations that only al- 
low us to observe the present day MF of this cluster. 
First, due to stellar evolution, the most massive star(s) 
have likely disappeared, and only stars with masses up 
to M maXj0 bs can be observed at the present day. Second, 
observational completeness (due to stellar crowding, for 
example) only permits observations of stars more massive 
than some lower mass limit, M comp . Third, the transla- 
tion of observed flux into stellar mass provides only an 
estimate of a star's true mass. That is, the 'observed' 
mass of a star could be an over- or underestimate of a 
star's true mass. Further, the observed mass of a star 
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Figure 2. A schematic representation of the Kropua IMF over the 
entire stellar mass spectrum for a hypothetical star cluster (Equa- 
tion[T|l. Here, jV* (red shading) represents the total number of stars 
in the cluster and N p[c ^ is the number of stars expected in the ob- 
served mass range bounded by M comp and M maXi0 t, s (blue shaded 
region), i.e., the MF. The upper end of a hypothetical power-law 
IMF is defined to be between M m i n and M max , which is the max- 
imal mass of a star that could have formed in the cluster. Mi; m 
is the upper bound of possible observed stellar masses, as set by 
stellar physics. The aim of this paper is to outline an approach 
for inferring the values of a, M max , and AT ple< j of a young cluster, 
given a set of N observed stellar masses. 

is bounded by M comp and Mn m , the maximum mass al- 
lowed by stellar evolution models, which may or may 
not reflect the upper stellar mass limit in nature. Each 
of these effects must be carefully considered in order to 
make an accurate measurement of the present day MF 
of a cluster, a necessary step to make any meaningful 
statements about a cluster's stellar IMF. 

For clarity in exposition, _we will begin with the sim- 
plest case where the data, d, consists of a set of N stars 
with perfectly known stellar masses, {Mi}, and the com- 
pleteness function is also perfectly known. Given the 
data, d, we aim to place constraints on the slope of the 
MF, a, the maximum mass of a star that could have 
formed, M max , and the expected number of observed 
stars, iVp re d, which normalizes the amplitude of the MF 
for this cluster (see Figure [2| . Throughout this paper, 
we assume that we have the observed completeness func- 
tion (either as a function of mass or luminosity) in hand. 
In practice, the completeness is a function of both the 
stellar fluxes (i.e., through the flux limit of the observa- 
tions) and the local image crowding, and is determined 
through extensive artificial star tests. Given that most 
cluster observations, certainly those beyond the Milky 
Way, are crowding limited, it will always be more chal- 
lenging to detect and accurately catalog a particular star 
if it is located near the center of a cluster rather than in 
the outskirts. For the purposes of this paper, we make 
the simplifying assumption that the completeness func- 
tion is independent of position. 

2.1. Deriving the General MF Posterior Probability 
Distribution Function 

We now lay out the general probabilistic framework 
for constraining the MF of a young cluster. We start by 
describing the MF of the cluster (cf. Equation [lj as: 

$(M) = N*p UF (M\6) , (3) 



where is the total number of stars in the cluster, 
Pmf(M I 0) is the probability distribution function (PDF) 
for the mass of any star in the cluster, and M is the true 
mass of the star; the MF is described by the parameters 
9 = {a, M max }, where a and M max refer to the high mass 
portion of the MF, i.e., a single sloped power law such 
as described by in Equation [Tl 

The PDF for the mass of an observed star, i.e., those 
that are massive (luminous) enough to end up in our 
dataset, is given by 

i \ PMF{M\e)p{obs\M) 

Pmf (M I 9, obs) = , (4) 

p(obs | 9) 

where p(obs | M) is the probability of observing a star 
given its true mass M, i.e., the completeness func- 
tion, and ^(obs | 9), the normalization necessary to make 
Pmf (M I 9, obs) a PDF. This normalization can be writ- 
ten as 

p(obs \9) = J pus{M I 9) p(obs | M) dM . (5) 

In this model, the expected number of stars in the data 
list, iVp r ed, sets the amplitude normalization of the MF, 
and can be expressed as 

/■oo 

iVpred = JV. / PMF a (M \ 9, obs) dM . (6) 

Jo 

However, given that we do not know a priori and 
that much of the MF is, as in many practical applica- 
tions, well below our detection limit, we cannot explic- 
itly solve for iVp re d, and therefore we will treat JV prc d as 
a model parameter that is independent of 9. 

So far, this model makes predictions for the probabil- 
ity distribution for the mass of any one observed star, 
Pmf (M I 9, obs), for a given set of parameters, 9, and 
for the expected number of observed stars, JV pre d- As- 
suming that each of the N observed stars in the cluster 
is drawn identically and independently from the under- 
lying mass distribution, we can write the probability of 
measuring the set of N masses, {M^}, as: 

N 

Pm ({Mi } I 9, {obs J) = n p MFo (Mi | 6, obs, ) . (7) 

1=1 

At the same time, we must also include the probability 
of actually observing N stars in the cluster, as N itself 
is a datum. This is the Poisson probability of observing 
N stars given an expectation of JV pre d stars and can be 
written as 

JV D ld e-^ d 

PPoisson (JV I JVpred) = • (8) 

We now have the two expressions that give the probar- 
bility of the data ({Mi}, N) given model parameters (9, 
JVpred) and the fact that data were observed, i.e., 'obs'. 
However, we wish to infer the inverse, namely the values 
of 9 and JV pre d given a set of N stellar mass measure- 
ments. Using Bayes's theorem, we can then write down 
the probabilities of the model parameters given the data, 
i.e., the posterior PDFs (pPDFs) for 9 and JV pre d as: 
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P P ost(0| {MJ,obs) =PM({Mi}\8,{obSi})p piioi (6) 



N 



= Y\_PMF ( M i I 9> °bSi) PpriorW 



(9) 



and 



Ppost (iVpred I N) Ppoisson (^l^pred) Pprior(-^pred) 
TN „-N pIcd 

Pprior 

(N pied ) , (10) 



N A . e 

pred 



where p pr ior(0) and Pprior(-'Vpred) reflect any prior or ex- 
traneous information or constraints on 9 and N pTed , re- 
spectively. In Equations [9] and [10| we have presumed 
that the 'evidence' terms, p pr ior({-™7}) and p pr i OT (N), are 
both constant, and have omitted them. In practice, the 
independence of these two pPDFs allows us to compute 
them separately, and then simply multiply them together 
to get the general pPDF for the MF: 



So far, we have assumed that the uncertainties on 
masses are negligible. If follows that a star's observed 
mass PDF is a delta function, and therefore it is directly 
interchangeable with a delta function PDF for the true 
mass of a star. Consequently, we can use Pmf d {M | 9, obs) 
as the generative model for the observations such that 



p MPo (M | 9, obs) = p(obs | M) cmfM M~ a = 

C M F o (0)M-< 

0, otherwise 



, M comp < M < M n 



(13) 



where the normalization cmf (9) is given by 



-MF, 



p(obs | M) M~ a dM 



Mr, 



M,., 



M- a dM . 
(14) 

In this form, the range of the MF model is set by the 
completeness limit on the low mass end, and the mass of 
the most massive star to have formed on the upper end. 

U sing the explicit expressions for each term in Equa- 
tion 11 we can now re- write the pPDF for this idealized 



case as: 



P P ost(0, N picd \{Mi},N, obs) = 
P pos t(9 | {Mi}, obs) Pp OS t(N pie d 



N). 



(11) 



The above derivation is limited to a single power- 
law slope, where the stars are well-above any break or 
turnover in the MF. Of course, in cases where data on 
stars well below 1 M Q are available, the same formalism 
carries through, just replacing a with a. 

3. DERIVING MASS FUNCTION CONSTRAINTS FROM 
IDEALIZED DATA 

We now show how to explicitly constrain MF param- 
eters given: (1) the pPDF from Equation 11 (2) evi^ 
dence for the data; (3) and any prior information on 
and N prod . In practice, the most informative prior on 9 
frequently comes from M max , where astrophysical con- 
straints may exist, especially if there is independent in- 
formation on the age of the population. Other general, 
although less physically informative, priors are require- 
ments that the number of stars is positive and that a 
does not take on unreasonably extreme values that have 
been ruled out by previous observations, e.g., a < —5. 

For the purposes of this initial idealized exercise, we 
adopt a simple 'boxcar' form for the probability of ob- 
serving a star, given its true mass (i.e., the completeness 
function) , 



p(obs | M) 



1, O.5M <M, 
0, otherwise 



JWeomp < M < Mii m 



(12) 



where Mn m is the maximum observable mass as set by 
stellar evolution, M comp is the minimum mass as deter- 
mined by the observational completeness limit. For sim- 
plicity, we will assume that the MF is a single power-law 
over the interval M comp < M < M lim and that M comp 
is the same value for all locations in the cluster; we dis- 
cuss the case of spatially varying completeness in S|6] In 
practice, p(obs | M) is the completeness function of the 
data. 



Ppost 

JV 

Ppoisson CMF o (0) M i a ) Pprior (9) Pprior (N pred ) , 

i = l 

(15) 

In practice we calculate the natural logarithm of the 
pPDF for computational ease. Also, note that in this 
idealized limit, p pos t{9, N plcd \{Mi}, N, obs) = in the 
event that M l < M comp or M max > M lim . 

3.1. Sampling the the pPDF with a Markov Chain 

To place constraints on the model parameters of inter- 
est, we sample the pPDF using a Markov chain Monte 
Carlo (MCMC) algorithm. MCMC techniques provide 
an efficient discrete sampling and clear interpretation 
of multi-dimensional spaces, such that the density of 
samples is highest around the most probable paramc- 
ters and lower in re gions of less probable values (e.g. 



Gelman et al. 19961. As a result, MCMC techniques 
are able to produce well-defined uncertainties for both 
the one dimensional (marginalized) distributions (e.g., 
p(a\ {M { }, N, obs), p(M max |{AfJ,A7,obs)) and reveal 
degeneracies between two dimensional (joint) distribu- 
tions (e.g., p(a,N prcd | {Mi}, AT, obs)). An MCMC ap- 
proach also permits efficient sampling in cases where 
more parameters are needed to adequately model the 
MF, e.g., such as when using a multi-component power- 
law to characterize data containing lower mass stars. 

In this analys is, we use the MCMC s a mpler emceq^"] 
as described in Foreman-Mackey et al. (2012). emcee 
is a pure-Python implementation ot the attine invari- 
ant MCMC en semble sampler developed by Goodman 
& Weare (2010). This sampler explores parameter space 
through a set of 'walkers'. For each MCMC increment, 



10 http://danfm.ca/emcee 
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each walker takes a step in parameter space by choos- 
ing another walker and moving along a line in parame- 
ter space that connects itself to the other walker. The 
size of the step is chosen stochastically, allowing for both 
interpolation and extrapolation, and the choice in step 
for each walker is based on the co-variance of the set of 
walkers. After each step, the pPDF is evaluated for the 
new set of parameters. Steps that increase the proba- 
bility are always accepted, whereas steps that result in 
a lower probability are sometimes accepted. For a suf- 
ficiently large number of steps, the ensemble of walkers 
samples parameter space with a frequency proportional 
to the pPDF. In a rough sense, this approach is akin 
to havi ng a set of parallel Met r opolis-Hastings MCMC 
chains (Metropolis et al. 1953 Hastings| 1970), but is 
much more efficient than Metropolis-Hastings sampling, 
as measured by the autocorrelation time (i.e., the time 
spent per function call; see discussions in |Goodman fc| 



|Weare|[20T0l |Foreman-Mackey et al.||20l2"l and the Ap- 
pendix of this paper) . 

emcee requires minimal initial input from the user. 
One must select a number of walkers and their initial con- 
ditions, and designate the number of burn-in steps and 
length of each chain. Followin g the recommendations in 



Foreman-Mackey et al. (2012) and experimentation with 
the simulated datasets, we chose 16 walkers and selected 
random (valid) values for each walker's starting point. 
After testing several combinations, we found that gen- 
erally 300 burn-in steps and 100 chain steps per walker 
resulted in stable solutions. To aide in computational 
efficiency, we placed sensible restrictions on the allowed 
ranges for each of the three parameters (i.e., p pT i 01 (9) and 
Pprior (^Vpred)): -6 < a < 6, M max , obs < M max < 120 
M Q , and < iV prc d < 10 2 x N. In the case of finite 
mass uncertainties, we r evis ed the restriction on M max 
for reasons discussed in |4.1| 

We refer the reader to the Appendix for a more in 
depth explanation of MCMC sampling theory. 

3.2. Illustrative Examples with Idealized Data 

We now illustrate this approach by applying it to ide- 
alized mock data (i.e., no mass uncertainties and a per- 
fectly known boxcar completeness function). A sim- 
ple way to create the simulated data for a power-law 
MF, assuming negligible mass uncertainties, is to de- 
fine a variable x via M — exp (ln[a; (1 — a)]/(l — a)) 
and then draw a uniform random variate in x between 

^min/max itinnt = (1 — a) _1 Af . . ,. For allilhlS- 
min/ max, input v j max/min, input 

trative ex amples in this paper we have adopted ainput ~~ 



2.35 (i.e., Salpeter||l955 ). We discuss the results of sim- 
ulations run with other values of a at the end of this 
section. 

Each dataset consists of N masses that have been dis- 
cretely sampled over the interval [M minMput , Af maXjinput ]. 
The value of M m i ni i nput is equivalent to the lower mass 
limit, Af comn , as specified by the completeness function. 



^comp 

For these exercises 



we have fixed AT, 



min, input 



= 3 M 



a conservative value f or resolved clusters in M31 (e.g. 



Dalcanton et al. 



2012). 



The mass of the most massive star that could have 
formed in a cluster is of great interest astrophysically. 
Essentially it provides a clue to the fundamental nature 
of how stars form, e.g., is there a universal value for the 
most massive that can form or does it depend on local 



environmental conditions? There are several studies in 
the literature that attempt to quantity the relationship 
between e nvironment and the m ass of the most massive 
star (e.g., Weidner et al.||2010 
dridge||20T2l r 



Cervino et al. 2011 El- 



In our simulations, the value of the most massive star 
that could have formed in a cluster is M maXj i nput . The 
mass of the most massive star "observed", i.e., on the 
mass list, is M maXiODS , which is less than M maXi i nput due 
to the effects of stellar evolution or dynamical ejection. 
Under the assumption that the mass of each star is in- 
dependently drawn from the IMF, it is fairly intuitive 
that M maX! obs cannot inform us about M maX! i nput beyond 
providing us a simple lower limit on Af maXimput , i.e., the 
most massive observed stars does not contain informa- 
tion about more massive stars that are no longer in the 
cluster. 

For a given dataset, we aim to recover the input MF 
slope, ainput, the maximal mass of a star that could have 
formed Ai maXimput , and require that our model produce 
a normalization that is consistent, i.e., within Poisson 
uncertainty, with the number of stars observed, N. We 
denote the corresponding model parameters as a, A/ max , 



and N, 



prcd • 



We applied this technique to mock data for different 
permutations of a m put, Af maXimput , and N. Results from 
the first example are shown in Figure 3} Here we show 
the recovered joint and marginalized distributions for a 
mock cluster which was generated assuming N = 1000 



stars and M„ 



Af„ 



= 54.0 



put = 60 M Q . We recover a = 2.35±g;g| 5 
j_ 8 , and JVpred = 995^31, where the indi- 
cated values represent the median and 0.5 x (84th - 16th 
percentile) of the corresponding marginalized distribu- 




Figure 3. The recovery of MF parameters a, M u 



N, 



pred in 



the case of idealized mock data, i.e., negligible mass uncertainties 
and a boxcar completeness function, from a single run of emcee. 
The joint distributions are shown as a scatter plot in the center of 
the first three panels, for the indicated parameters. Dark colored 
points are close to the most probable value, while light colors are 
farther away. The solid black line encloses 68% of the points cen- 
tered around the most probable combination of parameters. The 
corresponding marginalized distributions have been projected onto 
the axis of each plot. The solid red lines indicate the median of 
each marginalized distribution, and the shading encloses 68% of 
the distribution centered around the median. The roundness of 
each joint distribution indicates little degeneracy between the MF 
parameters. For this example dataset, we have recovered a = 
54.01ns, and 7V Drcd = 995±??. 



n oc + 0.05 M 



In the lower 



ative mass function is shown 
as the thick purple line, and 100 MF models randomly drawn from 
the pPDF are over-plotted as thin grey lines. 
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Figure 4. The recovery of the MF parameters a, M max , N pIc ^ 
in the case of idealized mock data, i.e., negligible mass uncer- 
tainties and a boxcar completeness function, from a single run of 

2.35, and 



*input 



emcee. This simulation has N = 500 stars, a u 
A^max, input = 60 . From top left to lower right: the marginal- 
ized distributions for a, M max , and 2Vp re d' The solid red line rep- 
resents the median of the distribution and the grey shaded region 
highlights 68% of the distribution centered around the median, i.e., 
the range enclosed by the 16th and 84th percentiles. In each panel, 
the blue dashed line represents the 'truth', i.e., aj nput , M max >0 bsi 
and N, respectively. For this example dataset, we have recovered 

_ 2 35"*" A<T — K1 q + „„J /\r . _ c; nn -(-22 + 



Mn 



■ 63.3_ 31 , and Af pro( j = 500_ 2 i> indicating 
excellent consistency with the input values. In the case of M max , 
we are only able to place constraints on M max |-> a and not on the 
mass of the most massive star formed. In the lower right hand 
panel, the observed cumulative mass function is shown as the thick 
purple line, and 100 MF models randomly drawn from the pPDF 
are over-plotted as thin grey lines. 



tion, i.e., l-cr for a Gaussian distribution. 

In general, the recovered values reflect the input val- 
ues with excellent accuracy and precision. However, as 
expected, the recovered value of M max is consistent with 
the mass of most massive observed star, M maXi0 b s , and 
not that of the most massive star that could have formed 
M maXj input- Thus, this constraint provides a lower bound 
on the upper mass limit of the IMF. 

The joint distributions in Figure [3] show little co- 
variance between parameters, which indicates that there 
is little degeneracy in this problem, i.e., the estimation of 
one parameter will not affect the accuracy and precision 
to which another parameter can be constrained. Conse- 
quently, and for clarity, we will focus on the marginalized 
distributions for the next set of examples. 

In Figures |4|6j we show select examples for tests of 
cluster MF recovery in the cases where N = 500, 100, 
and 25 stars and M maXi i nput = 60, 30, and 15 M Q . The 
first three panels in each figure show the marginalized 
distributions for a, M max , and iV pre( j. For the case of 
500 observed stars and Af ma x,obs = 59.9 M Q (Figure [4]), 
we recover a = 2.35+q'o7, -^max = 63.3l|"5, and -/V pred = 
500 j^, which are all consistent with their input or ob- 
served (in the case of M max ) values. 

Figures [5] and |6] arc examples of the same tests applied 
to populations with fewer stars and smaller dynamic 
ranges in mass. In both examples, the recovered parame- 
ters are consistent with the input values. As expected, as 
the amount of observational information decreases, the 
constraints become increasingly broad. However, even 
in the limit of 25 observed stars and M max i npu t = 15 
M (Figure [6], we find a = 2.50lg;^ and N pied = 26±t 




Figure 5. 

-^max.ijiput 

32. "T 



The same as Figure [4] only with N 
= 30 Mq. Here, we recover a = 2.37 



100 stars and 



+0.19 



prcd 



101 



-11 

-10- 



The recovered fractional precision on 



a is roughly a factor of 2.5 larger than for the example considered 
in Figure HI 



Thus, these ranges contain the input /observed values for 
each parameter. 

In general, these exercises demonstrate probabilistic 
approach in fj2] can recover MF parameters, even with 
limited information. In other words, useful constraints 
can be derived even for clusters that fail to reach a 'gold- 
standard' level. Importantly, we verify that the recovered 
MF slope is equivalent to the input IMF slop e, as is ex- 



pected in the case of a single age cluster (e.g., Elmegreen 
& Scalo 2006 ). This finding affirms that this probabilistic 
technique provides the means to directly infer the slope 
of the high mass IMF of a given cluster with full and ac- 
curate accounting of the associated uncertainties, in the 
absence of dynamical effects, unresolved binaries, etc., as 
discussed in £j6] 

On the other hand, the recovered constraints on Af max 
are less informative about the most massive star that 
could have formed, i.e., the upper mass limit of the IMF. 
We find that the recovered values of M max only constrain 
Af maXi0bs , which is only a lower limit on M roaXiillpu t , the 
IMF maximum mass limit for a given cluster. We there- 




is 20 25 30 35 40 45 50 



Stellar Mass [M : ] 



Figure 6. The same as Figure [4] only with N = 25 stars and 
Af m ax,mHut = 15 Mq. Here, we recover a = 2.50+°'gg, M max = 
13.6_ 2 3 Af prc( j = 26_ 5 . The recovered fractional precision on a 
is roughly a factor of 8 larger than for the example considered in 
Figure B] This exercise demonstrates that it is still possible to 
place infer the slope of the MF from a sparsely populated cluster. 
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-0.2 



-0.1 0.0 0.1 0.2 

{(^recovered) &input 



1.0 



1.5 



2.0 2.5 3.0 

Log N 



3.5 



Figure 7. The accuracy and precision to which we can recover a, in the case of negligible mass uncertainties and a boxcar completeness 
function. Left panel - The distribution of median values of a from 100 recovered datasets, each with 10 4 stars and M max j nput = 60 Mq. 
The resulting distribution is well-characterized by a normal distribution with fi = 0.006 and a = 0.07, indicating that this approach to 
MF recovery is free of systematic biases. Right panel - A 2D map of the theoretical precision to which a can be recovered as a function 
of the number of stars and dynamic mass range. This plot indicates that that intrinsic precision (i.e., 0.5x (84th-16th percentiles)) for the 
recovery of a ranges from < 0.1 in the limit of many stars and a high dynamic mass range to > 0.7 for a poorly populated cluster with 
a low dynamic mass range. However, under the assumption that masses are independently drawn from the MF, it is possible to combine 
multiple small clusters to improve precision of MF recover. The asymmetry in the distribution indicates a slight increase in the recovery 
of a for a small cluster with a high dynamic mass range relative to a well-populated cluster with a low dynamic mass range. 



fore suggest that dedicated searches for the most mas- 
sive star (e.g., spectroscopic surveys) are more suited 



than application of this technique 


; (e.g., Massey et al. 


1995a; Lennon 1997; 


Masscy 2003; 


Kudritzki et al.|:J0(W| 


Crowther et al.||2010 


Evans et al. 


201 1|). Consequently, 



the recovery of a, and only discuss constraints on M max 
in cases of particular interest. 

Extending the above exercises to larger ranges of N 
and M max allows us to illustrate the precision and ac- 
curacy to which a can be recovered for a diverse set of 
clusters. We first check the accuracy of the method to 
verify that there are no systematic biases in the recovered 
MF slope. To do this, we simulated 100 different clusters 
each with 10 4 stars, to minimize stochastic effects, and 
A^max, input = 60 M , for the same underlying MF. We 
recovered the MF for each of the 100 clusters and plot 
the distribution of the difference between ai npu t and the 
median value of ^recovered in the left panel of Figure [7J 
The mean of the distribution is 0.006 and 0.5 x 84th - 
16th percentiles, i.e., a Gaussian 1-er, is ±0.07. Based 
on this set of tests, we see no evidence for systematic 
biases in the recovery of a in any of the parameter space 
considered. 

Next, we explore the theoretical precision to which a 
is recovered as a function of N and the dynamic range of 
the observed masses, as proxied by log(M maXi0 b s /Af comp ). 
We constructed a grid of model clusters spanning a range 
in N and log(Af maXj0 bs/Af comp ) as follows. We fixed 



ture of the sampling, we applied a Gaussian smoothing 
kernel to the resulting grid with a l-cr width of ~ 5% of 
the dynamic range in each variable. We computed the 
smoothed and un-smoothed values and found they differ 
by no more than ~ 7% at any point. 

We have plotted the results of this exercise in the right 
panel of Figure [7j In the limit of high observational in- 
formation, e.g., N ~ 10 3 stars and log(M maX;0bs /M comp ) 
~ 1.1, we see that a is recovered with a precision of 
~ ±0.1. At the other extreme (N ~ 10 stars and 
log(M max , obs /Mcomp) ~ 0.3), the precision in a roc0 vcrod 
decreases significantly to > ± 0.7. Because the mock 
data used for these tests are highly idealized, i.e., per- 
fectly known individual stellar masses and completeness, 
these values reflect the highest theoretical precision with 
which a can be inferred from a given number of stars and 
a measured dynamic mass range. As we discuss in the 
next sections, other considerations such as finite mass 
uncertainties and completeness will degrade these levels 
of precision, lessening the achievable constraints. 

Additionally, using a least squares fitting routine, we 
find that Aa can be reasonably approximated as a func- 
tion of log(iV) and log(M maXj0bs /M comp ) by a 3rd order 
polynomial 



A« = ^]T/Mlog(A0y Qog(M, 



max.obs 



/Mi 



comp 



)) j 

(16 



comp 



M, 

drawing random values N and Mmax.input such that 0.8 
< k>g(A0 < 3.5 and 0.2 < log(M maX; obs) < 1-3, and dis- 
cretely sampled the IMFs following the procedure out- 
lined at the beginning of this section. We recovered the 
MF slope for each realization and computed Aa = 0.5 x 
(84th - 16th percentile), i.e., a percentile based l-cr equiv- 
alent. To account for fluctuations due to the discrete na- 



= 3 M Q and created 7000 simulated clusters by where the coefficients ftj are listed in Table 

Within the range of log (AT) = 1.5 to 3.0 and 
log(M maXi0b s/M comp ) = 0.4 to 1.2, this approximation 
is generally good to within ~ 20%. Outside this range, 
this analytic expression is noticeably steeper, meaning 
the extrapolations into regimes of lower and higher in- 
formation content lead to over- and underestimates of 
the true precision. Nevertheless, when used with a rea- 
sonable level of caution, this analytic expression provides 
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Figure 8. Left -The 2D map of the theoretical precision to which a can be recovered from Figure]^ We have over plotted the indicated 
precision from the literature studies shown in Figure [l] and tabulated in Table [T] Points that appear darker than the surrounding regions 
have listed errors above the theoretical precision limits, where as lighter points have errors that are smaller than are theoretically permissible. 
Right - A histogram showing the distribution of the ratio of quoted literature errors and our computed theoretical limits for clusters show 
in Figure [8] and tabulated in Table ^ The red shaded region indicates the region in which literature studies have quoted errors less than 
the theoretical limit. The range of ratios varies over an order of magnitude, and approximately ~ 3/4 of the literature studies quote error 
bars than are below theoretical expectation. 



a good rule of thumb for the best attainable precision on 
the MF slope. 

Finally, we tested the accuracy and precision to which 
other a values can be recovered. Differing values of a 
dictate the relative distribution of stellar masses, and a 
significantly steeper MF slope would produce fewer mas- 
sive stars, potentially resulting in less accurate and/or 
precise constraints on the MF slope. To test this po- 
tential dependence on a, we selected two different, but 
plausible values of a = 1.8 and 2.8, and repeated all the 
simulations presented in Figure m We found that both 
the accuracy and precision to which a can be recovered 
is statistically consistent in all cases. There were random 
variations at the few percent level; however, such small 
differences are easily accounted for by the discrete nature 
of the sampling and the finite number of realizations. 

3.3. Comparing Theoretical Precision Limits with 
Values from the Literature 

As an informative exercise, in the left panel of Figure [8] 
we have over-plotted the error bars of the literature MF 
values from Figure [l] on top of our derived theoretically 
attainable values. The contrast in colors between the 
points and shaded contours indicates the level of agree- 
ment between the two; points that are lighter than the 
surrounding region are below the theoretically attainable 
precision, while darker points are larger than this lower 
limit. A cursory inspection of this figure indicates a full 
range of differences between the two datasets, i.e., there 
are a number of points with extreme light and dark color 
contrast. We have tabulated the literature data used to 
construct this plot in Table [T] 

To better quantify differences in the two datasets, in 
the right panel of Figure[8]we have plotted a histogram of 
the ratios between the theoretical and literature MF pre- 
cision values. The histogram shows that the ratios vary 
over an order of magnitude around unity, with ~ 3/4 
of the error bars quoted in the selected literature stud- 
ies being smaller than the theoretically attainable limit 
(red-shaded region). The underestimated error bars are 



typically lower than the theoretical expectations by a fac- 
tor of ~ 2 or more. Overly small error estimates can give 
the false impression of a good constraint, which could 
lead to claims of a significantly unusual MF measure- 
ment. Conversely, large overestimates of the errors on 
a MF slope indicate that maximal information may not 
being extracted from the data. Both cases reinforce the 
importance of developing a well-defined methodology for 
correctly characterizing the MF slope of a young cluster. 

Using the literature measurements of a in tandem with 
the theoretical lower limit errors, we compute a mean 
a for the literature values considered. In Figure [9j we 




Figure 9. A modified version of Figure ^ based on a comparison 
between the amplitude of the error bars quoted in the literature and 
the theoretical lower limits derived in this paper. The blue points 
are studies whose quoted errors are larger than the theoretical lower 
limits. The magenta points are those that had quoted error bars 
smaller than the theoretical lower limit. For these points, we have 
assigned and plotted new error bars equal to the theoretical lower 
limits. We then use all the literature data to compute the weighted 
mean (solid green line) and weighted 1— a (lightly shaded green 
region) of a, which we find to be (a) = 2.46 with a l-cr dispersion 
of 0.35 dex. Note that while we are able to assign more robust 
error bars to most points, we are unable to assess the accuracy of 
each value of a, which would require re-analysis of each dataset. 
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have re-created Figure [T] but have plotted the theoretical 
lower limits errors for those studies which have underes- 
timated the uncertainties. These are indicated as ma- 
genta points. Studies whose error bars are larger than 
the theoretical limit have been left as is, and are plot- 
ted as blue points. From the Figure [9] and Table [T] 
we see there is no discernible trend for which studies, 
values of a, or mass ranges have underestimated error 
bars. Using all reported a values, along with the new 
set of uncertainties, we computed a weighted mean and 
weighted 1-er uncertainty for this selection of IMF stud- 
ies and find (a) = 2.46 with a \-a dispersion of 0.35 dex. 
This value is consistent with the slope of a Salpeter and 
Kroupa IMFs, a = 2.35 and 2.3, respectively. However, 
the scatter is large ensuring the mean value is als o con- 
sistent with the Kennic utt IMF (Kennic utt 1983), a = 
2.5, and the Scalo IMF ( |ScaIo||l986| ), a ~ 2.6. Although 
our computed mean value ot a depends on our selection 
of literature studies and the accuracy of each study, the 
broad scatter reinforces the notion that our knowledge of 
the high mass IMF slope is far from secure. 

4. INCORPORATING UNCERTAINTIES IN STELLAR 
MASSES 

Up to this point, we have explored placing constraints 
on MF parameters presuming data with negligible in- 
dividual mass uncertainties. However, in practice, the 
mass of a star is not perfectly known. Degeneracies with 
extinction, metallicity, etc., along with uncertainties in 
stellar models all contribute to uncertainties in estimat- 
ing the mass of a star. 

Clearly, the uncertainty on the mass of an a star needs 
to be incorporated into the pPDF for realistic MF deter- 
mination. To understand the relationship between MF 
determination and stellar mass uncertainties, it is useful 
to review how stellar masses are measured. 

Most individual star masses are inferred from photo- 
metric spectral energy distribution (SED) fitting, i.e., by 
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Figure 10. A schematic demonstration of the importance of con- 
sidering the MF as a prior in the determination of individual stellar 
masses.^The result of photometric SED fitting for a stellar masses, 
i.e., p(d\ M), is denoted by the blue line. The MF with a = 2.35 
(red dot-dashed line), indicates that lower stellar masses are more 
likely than higher stellar masses. Therefore, the best mass esti- 
mate for the star comes from the convolution of these two func- 
tions (green dashed line). Not taking the MF into account when 
deriving stellar mass PDFs can lead to noticeable overestimates of 
the star's true mass. 



comparing photometric observations of a star, d, with 
predicted fluxes from stellar models. Specifically, one 
can construct a set of predicted fluxes in the observed 
bands by for various combinations of log(T e //), \og(g), 
Z, Ay, M, age, etc., and then perform a comparison, 
e.g., ^-minimization, between the observed and pre- 
dicted fluxes. From this process one is left with a multi- 
dimensional distribution for the probability of the star's 
parameters given the data, i.e., p(di | log(T e //), \og(g), 
Z, Ay, M, age). One can then compute the likelihood 
of the data given a mass, i.e., p(d \ M), by marginalizing 
over all other model parameters. The resulting PDF for 
a star's mass accounts for correlations and uncertainties 
in the all other parameters. 

While the quantity p(d | M) is often thought to reflect 
the probability of a star's mass, the desired quantity is in 
fact, p(M | d), the probability that a star has a particular 
mass given the observations. To determine this quantity, 
one can apply Bayes's theorem, which yields 

p(M | d) = l -p(d\ M)p prioT (M) , (17) 

where j? pr i O r(-^0 reflects any prior knowledge about the 
probability of individual stellar masses, and C is a nor- 
malization factor needed to make p{M \ d) a true proba- 
bility. 

A temptingly simple assumption for p pr i O r(-^0 would 
be that all stellar masses are equally likely, ^ i.e., 
Pprior(-^0 = constant. In this case, p(M \ d) and p(d \ M) 
can be treated interchangeably, which is at least implic- 
itly assumed i n many published analyses (e.g. 



|Massey| 

et al.|1995a|b [ |Panagia et al.|[2000t [Brandner et aLPtiSj 
Bianchi et al. 12012}. 



However, it is clearly more appropriate to adopt the 
cluster's MF as a prior on the distribution of individual 
stellar masses. As illustrated in Figure [TUJ failure to 
consider the effects of the MF when determining the mass 
of a star will lead to a biased mass relative to the star's 
true mass in the presence of significant mass errors. 

Adopting the MF from Equation [4] as a prior on stel- 
lar mass, we can now write down an expression for the 
probability of the data, d, given MF parameter, 9, as: 



p(d\ 9, obs) = J p(d\ M) p MFo (M | 0, obs) AM 

_ p{ohs\d) J p(d z \M)p M F(M | 6) dM 



and 



p(obs | 6>) 



p(obs \0)= p(obs | d) p(d \6)dd, 



, (18) 



(19) 



where p{d \ M) is the likelihood of the observed data 
given a true masSj as determined by photometric SED fit- 
ting, andp(obs | 9) is the normalization factor, which rep- 
resents the probability of observing the data (i.e., fluxes) 
given the MF model parameters 9. 

Using Bayes's theorem, we can write the probability for^ 
the parameters, 9 and -/V prec i, given the observed data, d 
as: 
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likely to be constant. Stars with M > 25 M Q are less 



Ppost (ff, Npred I d, obs) = p poi (N | JV pred ) X 
^ p(obs \dj)J p(d t | M) pmf (M | 0) AM , 

/ . -t Pprior \y ) /'prior \i * prcd / j 

i=i P(°bs | 6) 

(20) 

which is the pPDF in the case of non-negligible mass 
uncertainties. 

4.1. Illustrating the Effects of Log- Normal Mass 
Uncertainties 

As an example of how to fold in finite mass uncertain- 
ties in practice, we consider the case that each star has 
an observed mass PDF that can be represented by a log- 
normal function. The assumption of a log-normal not 
only provides an analytically tractable solution, which 
aides in the clarity of this paper, but also log-normal 
approximations (i.e., log(M) ± 5m) have been used in 
the literature to summarize the inferred masses of stars 



e.g. JMassey fc Hunter|1998||Massey|2011[|Bianchi et al 



20121 . 

Consider a set of N stars each of which has an asso- 
ciated mass PDF. We believe that each observed mass 
PDF is the convolution of a single true mass, Mj, with 
a log-normal noise model (due to observational uncer- 
tainties, the flux-to-mass conversion process, etc.). The 
inferred mass for each star can then be summarized as 
having a mean value of Mj and a fractional mass uncer- 
tainty, fi = Mi /Mi, i.e., the linear width of the mass 
PDF scales as ~ e* . Considering only the upper portion 
of the IMF, the stars we observe can have true masses 
ranging from M m - m to M max , while the inferred masses 
will range according to the stellar evolution models and 
the observed completeness function, and thus have lim- 
its of Ar comp and Mn m . For simplicity, we will adopt a 
boxcar completeness function written in terms of M, 



p(obs | M) 



1, M comp < M < Mum M Q 
0, otherwise . 



(21) 



Having defined the relevant variables and the ranges, 
we can return to the derivation. The goal of this pro- 
cess is to place constraints on the MF parameters, given 
the set of log-normal PDFs for the inferred masses, 
and Equation 18 provides the framework to complete 



this task. Substituting a log-normal mass model for 
p(di | M) and making the coordinate transformation X 
I n M/M, allows us to write: 



p(di 1 0, obs) = 

(22) 

where _the completeness is over the observed masses, 
p(obs | M), the integration is over the permissible range 
of true masses, and we have made the simplifying as- 
sumption of a constant fractional mass uncertainty for 
all masses, i.e., fi — > /. This assumption has been made 
for clarity in this derivation, but in practice, /, is not 



Massey 



constrained that those of lower mass stars (e.g 
et al.|1995a ), owing to both uncertainties in massive star 
evolution models, as well as degeneracies in the optical 
and near-UV photometric colors of massive stars, which 
make it difficult to precisely characterize extremely mas- 
sive stars. 

Integration of Equation [22] yields a closed form solu- 
tion: 



p(d\6, obs) = 

c D (0) c MFo (0) M~ a p(obs j M) e 



(a-l)/ 2 +ln(^p) 



/ 



where 

cntf)- 1 = 
cmf o (0) e 

(a-l)/ a +ln0 w 



( a -i)f + H^) 



(23) 



M lin 
f2 i i„/M m 



M a p(obs|M)x 



f 



f 



dM , 



(24) 



and X] is the normalized cumulative distribution function 
of a standard Gaussian PDF with fi = and a = 1: 



V (X) 




dx 



(25) 



In the limit of very accurate mass estimates, i.e., / — > 
0, the expression for p(di\9) reduces to a scenario where 
the masses are perfectly known, as in Equation |13| 

Having derived an expression for p(di \9), we can now 
write the pPDF for the MF parameter as: 



Ppost (0, Nprcd I {di }, N, obs) 

Y 

Ppoi (N\N pled ) ^lnp(^|0,obs) Pprior (0) Pprior (-/Vpred) , 
1=1 

(26) 

4.2. Illustrative Examples with Mass Uncertainties 

In this section, we present several examples of MF re- 
covery in the presence of non-negligible mass uncertain- 
ties. However, before delving into specific cases, it is 
instructive to first examine how log-normal mass PDFs 
affect the observed MF. 

4.2.1. Schematic Examples with Continuous Mass Functions 
In Figure 



11 



we illustrate the changes in the (continu- 
ous) observed MF for different values of /. For reference, 
the solid black line represents a Salpeter IMF from 3 to 
120 Mq. The other plotted lines represent MFs gener- 



ated using Equation 22 for select values of / (0.1, 0.5, 
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Figure 11. A schematic demonstrating the influence of log- 
normal individual stellar mass PDFs on the MF. Each distribution 
was constructed assuming M maXj i npu t = 40 Mq (grey line) and 
"input = 2.35 (black line). We then varied the value the fractional 
mass uncertainties, /, as indicated. Small fractional mass uncer- 
tainties lead to a slightly steeper than Salpeter MF slope, and pre 
dieted the presence of masses larger than Afmax input- 
values of /, e.g., 1, the observed MF slope is flatter than Salpeter, 
with observed values of M max significantly larger than M max i npu t. 
Failure to properly account for mass uncertainties when modeling 
the obser ved MF can lead to biases in the recovery of a and M max . 
See 94.1|for more discussion. 



1.0, 2.0). In each case, Af maXi i n p U t = 40 M and each 
MF was normalized to have the same number of 3 M 
stars as a Salpeter MF. 

For the smallest fractional mass uncertainty, / = 0.1, 
we see that the observed MF (solid blue line) is identical 
to Salpeter until ~ 37 M . At this point it begins to 
deviate from the Salpeter MF, and predicts the presence 



of stars above the nominal value of M, 



max, input 



40 M^ 



Taken at face value, this MF appears slightly steeper and 
it implies a finite probability that M maX! ob s is larger than 
-Mmax, input- The MF for f — 0.5 exhibits the same qual- 
itative behavior, with an even steeper slope and a signif- 
icantly larger apparent value of M maXj0 bs- Naive fitting 
of the observed MF, without taking mass errors into ac- 
count, would therefore lead to biased measurements of 
both a and -M max . 

For the extreme values of f=l and 2, the observed 
MFs flatten out relative to Salpeter, and each predicts 
a significant population of stars with masses in excess of 
100 Mq. For such extreme uncertainties, a significant 
number of stars with true masses below 3 M will have 
inferred masses above the completeness limit. 

This schematic also illustrates the need to alter the 
prior restrictions placed on M max . In the case of no 
mass uncertainties, the inferred and true masses are iden- 
tical, implying that M max > M maX)0 b s was a reason- 
able requirement. However, for finite mass uncertainties, 
Mnax.obs is likely an overestimate of M maXtinput , i.e., the 
maximal inferred mass is larger than the maximal true 
mass, invalidating the previous assumption. For the ex- 
ercises involving finite mass uncertainties, we therefore 
only require that M max < M\\ m , the larger mass permit- 
ted by stellar evolution models. 

4.2.2. Illustrative Examples with Simulated Clusters 

We now use mock data to explore recovery of the MF 
in the case that each mass has a log-normal PDF. To 
generate the simulated data, we construct a convolution 



between the MF a log-normal noise model with a — 2.35 
(Equation 22 1. We then draw 10 6 stars from this func- 
tion, apply the boxcar completeness function from Equa- 
tion [21] and randomly subsample the mass list to obtain 
the desired number of stars. 

In Figures [12] and [13] we show results from the re- 
covery of a for select simulated clusters. In each case, 
the mock data consisted of 10 4 stars (to minimize ran- 
dom noise effects) and discrete values of / = 0.1, 0.5, 
and 1.0. To illustrate the effects of the observed dy- 
namic mass range, we also considered two different val- 
ues of M maXiinput = 20 M and 60 M . As a baseline 
for comparison, we applied two versions of the code to 
each dataset, one that included modeling of the mass un- 
certainties (the 'noisy' MF model; Equation 26 ) and the 
simple model MF (the 'noiseless' MF model; Equation 
|15[ ). In each panel, the dashed line indicates the input 
value of a while the solid line indicates the median value 
for each marginalized PDF. 
In Figure |12[ we show the recovery of a for 
For Extreme Afmax. input = 60 M . In the case of small mass errors (/ 
= 0.1), both the noisy (red) and noiseless (grey) model 
MF models provide excellent recovery of a m p Ut . The 
small level of offset between the results are consistent 
with the expected scatter (~ ± 0.03) in the case of N = 
10 4 stars. The two methods also return the same level 
of precision on a reC overedi as indicated by the comparable 
widths of the two distributions. 
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Figure 12. The recovered marginalized distributions for a when 
accounting (colored histograms) and not accounting (grey his- 
tograms) for individual stellar mass uncertainties in the MF model 
(colored histograms) and not including mass uncertainties (grey 
histograms). Each star has a log-normal mass PDF with a frac- 
tional uncertainty, / = M/M, whose assumed constant value is 
indicated. For each value of / the noisy (i.e., accounting for mass 
uncertainties) and noiseless (i.e., not accounting for mass uncer- 
tainties) models were applied to the same data, which always had 



max, input 



60 MQand N = 10 4 stars. In each of the first three 
panels, the true MF slope (Salpeter), is indicate by the dashed 
black line. Application of the noisy model yields recovery of the 
colored histograms, whose median value is indicated by the solid 
colored line. The grey histograms are the result of applying the 
noiseless models, and the solid black line represents the median of 
this distribution. For reasonable fractional mass uncertainties, e.g., 
/ < 0.5, there are only small differences in the MF slopes recovered 
by the two methods. For the extreme case, / = 1.0, modeling the 
mass uncertainties yields near perfect recovery, while failure to do 
so results in a systematically flatter MF slope. In the lower right 
panel, we show the MF distributions used in each case, and have 
included a Salpeter MF (with an arbitrary vertical scaling applied) 
for reference. 
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Figure 13. The same as Figure 1121 only with Mj 



1 max, input 



20 



Mq. In this case, the differences in the recovered MF slopes be- 
tween the two methods are more pronounced, even for modest val- 
ues of /. These examples demonstrate the importance of including 
a good characterization of mass uncertainties in order to measure 
the MF slope. 



For significantly larger mass errors (/ = 0.5), both the 
noisy (green) and noiseless (grey) models still provide 
excellent recovery of a. The distributions again share a 
comparable width. Although the noiseless PDF appears 
shifted toward slightly higher values of a, the median re- 
covered value is within the expected scatter for a cluster 
with N = 10 4 stars. 

Extending this exercise to extremely large mass errors 
(/ = 1.0) reveals significant discrepancies between results 
from the noisy (blue) and noiseless (grey) models. Based 
on the schematic in FigurcfTT] the observed MF should be 
flatter than Salpeter, which is what is the noiseless model 
returns, with the median a rcC ovcrcd = 2.20. Although the 
distribution is comparable in width to those in the other 
panels of the plot, it does not overlap the true value of 
a. In contrast, application of the noisy model recovers 
a to excellent accuracy, although to a lower degree of 
precision. 

To test the sensitivity of MF slope recovery to M max , 
we conducted the same exercise, only with Af max , input = 
20 Mq . T he recovered distributions for a are shown in 
Figure 13 For the lower value of M ma x,input> we see more 



pronounced differences in the results from the noisy and 
noiseless models. For both / = 0.1 and 0.5, the noiseless 
model recovers systematically higher values of a rela- 
tive to the input. While only a modest bias is present 
in the case of / = 0.1, the offset is quite substantial 
when / = 0.5. In contrast, the recovered values of a 
from application of the noisy MF model are in excellent 
agreement with a m p U t. Once again, the widths of the 
noisy model distributions are slightly broader than the 
noiseless model results. However, the slight decrease in 
precision is a small tradeoff relative to inferring a signif- 
icantly more accurate MF slope. 

In the case of / = 1.0, application of the noiseless MF 
model to the data yields an only slightly biased (< 0.1) 
median value of a. This is perhaps somewhat surpris- 
ing, given the extremely large fractional uncertainties. 
However, this small amount of disparity is purely coinci- 
dental. The combination of / and M maXj i nput conspired 
to produce an observed MF slope that appears close to 
the true MF slope. Such a case might lead to the conclu- 



sion of a precise constraint on a MF that is statistically 
different from Salpeter, despite the large mass errors; 
however, such a conclusion would clearly be false. In 
contrast, modeling the mass uncertainties results in an 
accurate recovery of the true MF, albeit with decreased 
precision, as previously discussed. 

Overall, we find that incorporating the effects of mass 
uncertainties into MF modeling results in an unbiased 
recovery of the MF. In contrast, failure to model even 
moderately large mass uncertainties, e.g., / > 0.5, can 
lead to biases that are difficult to quantify, which can un- 
dermine a physically meaningful interpretation of the re- 
sults. Additionally, neglecting to model mass uncertain- 
ties can lead to a severely low estimate of M max , which is 
equivalent of underestimating the cluster's age. Finally, 
as previously mentioned, in practice is unlikely to be 
constant, as the inferred properties of stars with M > 25 
Mq are far more un certain than lower mass stars (e.g., 
Massey et al.|[1995a[ ). Including a non-constant value of 
/ is straightforward in the presented framework. 

5. INCORPORATING THE EFFECTS OF OBSERVATIONAL 
COMPLETENESS 

5.1. Preface 

We now consider the effects of observational complete- 
ness on the inference of MF parameters. For clarity in 
the mathematics and the examples, we return to the case 
of perfectly known masses. In practice, both complete- 
ness and mass uncertainties need to be simultaneously 
modeled, but they typically affect opposite ends of the 
observed mass spectrum, meaning that they are largely 
not degenerate in terms of the data they impact. 

When deriving the general pPDF for a cluster's MF 
(Equation |11[ ) , we implicitly included the completeness 
function ihthe probability of a star being on our data list, 
Pmf (M I 0, obs). However, writing the pPDF to explic- 
itly include the effects of completeness provides a better 
illustration of the role it plays. Using Equation[4j we can 
write the pPDF for a general completeness function as: 



P P o S t(6», N pled I {Mi}, N, obs) = p poi (Af | N plcd )x 

N 

J2(p(° bs I M ) C MF o 0) M~ a ) PpriorWPpriorCAWd) . (27) 
i=l 

Here we see that the completeness enters into the pPDF 
as a linear weight on the MF. In the limit that the 
p(obs | M) — constant, it will have no effect on the pPDF 
as all points are equally weighted. In the next section we 
illustrate how a more realistic completeness function af- 
fects the recovery of the MF. 

5.2. Illustrating the Effects of a Linear Completeness 
Function 

In practice, a boxcar function is not a realistic rep- 
resentation of observational completeness. Complete- 
ness functions in clusters vary widely in form depending 
largely on the surface brightness of the cluster and the 
resolution of the observations. In this section, we con- 
sider a linear ramp completeness function. Such a lin- 
ear completeness function is analytically tractable and 
provides a reasonable first order approximation to more 
realistic c ompleteness functions presented in the litera - 
ture (e.g., Anderson et al. 2008 Gennaro et al. 2011). 
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The linear ramp function, which ranges from to 100% 
completeness, can be written piecewise as 



p(obs | M) 




M > Mcmax 

M < M CmiD 

-^Cmin < M < Mcmax , 

(28) 

where Mc max and -Mcmin are the upper and lower bounds 
of the linear portion of the completeness function. 

We again turn to mock data with no mass uncertainties 
to illustrate MF recovery in the case of a linear complete- 
ness function. Specifically, we consider the cases of how 
perfectly known, slightly incorrect, and extremely incor- 
rect knowledge of the completeness function affects MF 
recovery. 

To genera te th e mock data we followed the general de- 
scription in §3.2[ with a few modifications. After drawing 
a list of masses from the power-law MF, we applied the 
'true' completeness function to the mass list. Motivated 
by data from PHAT, for this exercise, the true complete- 
ness function is always a linear ramp between Mc m m = 
2 and Mc max = 4 M Q , with a completeness fraction of 
for stars with M < 2 M and 1 for M > 4 M Q . Af- 
ter applying this completeness function to the mass list, 
we then rejected stars from the list with a probability 
equal to their completeness fraction, e.g., a mass with 
j»(obs | M) — 0.3 has a 30% chance of being observed. We 
then perform a uniform draw from the remaining stars to 
arrive at the desired number of stars. This set of masses 
constitutes the true mass list. To examine effects of a 
non-perfectly known completeness function, we followed 
the same procedure as descried above, only the values 
of Mcmin and Mcmax are different when applied to the 
mock data and when used to model the completeness in 
the MF recovery. 

Keeping the same true completeness function, we then 
explored MF recove ry fo r various assumed completeness 
functions. In Figure 14 we considered the case in which 
the true (grey) and assumed (red) completeness functions 
were identical linear ramps between 2 and 4 M Q . We sim- 
ulated 100 datasets each with 10 5 stars and M maX! i nput = 
60 Mq, and measured the median value of a re covered) 
whose distribution is plotted as the red histogram in the 
right panel. The median of this distribution is repre- 
sented by the solid red line. For a perfectly known com- 
pleteness function, we found excellent recovery of the MF 
slope, as expected. 

We next considered the cases where the true and as- 
sumed completeness functions were not identical. In the 
first case, the ramp in the assumed completeness func- 
tion occurred between Mc m in = 2 and Mc max = 3 M 
(green), meaning the data are assumed to be more com- 
plete actually are, i.e., the stars between 2 and 4 M Q are 
assigned artificially high weights in Equatio n [27} As be- 
fore, we simulated 100 different clusters of 10 stars, and 
examined the distribution of the median recovere d va lues 
of a. As shown in the right-hand panel of Figure [l4j the 
recovered values of MF slope are clustered around a = 
2.0, with small scatter. This test indicates that a moder- 
ate overestimation of the true completeness function can 
induce a significant underestimate of a, because stars 
that are "missing" from the lower mass end due to ob- 
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Figure 14. The recovery of a from from various linear complete- 
ness functions. Left Panel - Modeled linear completeness functions 
that are perfect (red dotted), moderately underestimated (blue 
dashed), and moderately overestimated, compared to the "true" 
completeness of the data (grey). In this example the true com- 
pleteness is a linear from from 2 to 4 Mq, with a value of 0% 
below 2 Mq and 100% above 4 Mq. Center Panel - A histogram 
of the data (grey) with the "true" completeness function applied. 
The case of perfectly modeled completeness is in red. A moder- 
ate underestimate of the completeness (blue) indicates that more 
lower stars are present than expected, whereas a moderate overes- 
timate of the true completeness (green) indicates a deficit of lower 
mass stars relative to observations. Right Panel - A series of his- 
tograms indicating the recovered median values of a for each of 
the modeled completeness functions. The distribution of median 
values is derived from 100 different datasets, each with 10 stars 
and M max i nput = 60 Mq. The solid colored lines indicate the 
median of each distribution, and grey-dashed line is the di npu t. 
A well-characterized completeness function (red) results in an ac- 
curate recovery of a, whereas a modest overestimate (green) or 
underestimate (blue) lead to systematic under and overestimates 
of a. 



servational completeness are instead interpreted as being 
underp roduced by a flatter M F. This finding is similar to 
that in |Ascenso et al. ( 2009[ ), who demonstrates that an 
artificial flattening of the MF due to completeness con- 
siderations mimics the behavior of predicted mass segre- 
gation effects. We discuss this point further in Sj6] 

We also considered a case of a modest underestimation 
of the true completeness function. Here, the mass lim- 
its on the assumed completeness function were Men 
2 Mq and M Cmax 



6 Mq. In this scenario, stars be- 
tween 2 and 6 Mq are assigned artificially low weights in 
Equation [27] Application of this assumed completeness 
function results in median recovered values of a = 2.7 



(blue in Figure 14). This exercise indicates that even a 



moderate underestimation of the true completeness func- 
tion can induce a significant overestimate of a. 

Some MF studies in the literature cither apply a con- 
servative completeness limit, i.e., only consider data that 
are likely to be 100% complete, or make no attempt to 
correct for completeness. We explore the effects of each 
of these scenarios and present the results in Figure \TE\ 
In the case that no attempt is made to account for com- 
pleteness (green), i.e., all data above 3 M are complete, 
the recovered slope of the MF is ~ 1.7, which is signifi- 
cantly flatter than the true MF slope of 2.35. 

The scenario in which a conservative compl eten ess cut 
is applied is more promising (blue in Figure 15). Here, 
we see that the input value of a is accurately recovered, 
albeit with lower precision as both the number of stars 
and dynamic mass range are smaller. 

In this section, we have shown that by correctly char- 
acterizing completeness, through an example with a lin- 
ear completeness function, the MF slope can be recov- 
ered without bias. However, failure to include well- 
characterized completeness corrections can lead to sys- 
tematic biases. This finding can readily be generalized to 
more complex completeness functions, which simply in- 
volves replacing p(obs | di) by another functional or tab- 
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Figure 15. The same as Figure [H] only for extreme completeness 
examples, i.e., all data are incorrectly assumed to be complete 
(green) and a conservative completeness cut (blue). The "true" 
completeness function (grey) is a linear ramp from 2 to 4 Mq. In 
this case, the assumption that all the data are complete leads to 
a large deficit of lower mass stars (green, center panel), resulting 
in an severe systematic underestimate of a (green right panel). 
A lack of appropriate treatment of completeness has significant 
implications for und erstanding mass segregation in cluster (e.g., 
|Ascenso et al.|2009| . In contrast, a conservative completeness cut 
(i.e., a step function at 8 Mq ) , permits an unbiased recovery of the 
MF slope, at the expense of increased scatter due to the smaller 
dynamic mass range and decreased number of stars. 



ulated form. To determine a completeness function in 
practice, we strongly advise that extensive artificial star 
tests should be performed. Although this process can 
be a computationally expensive, it is essential to have a 
well-characterized completeness function in order to ac- 
curately measure the MF slope. 

6. TOWARD FULL FORWARD MODELING OF A 
CLUSTER'S MF 

Throughout this paper, we have employed a simple 
model of the IMF, and not considered other physical 
(e.g., dynamics) and observational effects (e.g., unre- 
solved binaries), which can influence the accurate mea- 
surement of a MF from resolved stars. Here, we briefly 
discuss how such outstanding issues can potentially be 
incorporated into the context of a probabilistic frame- 
work for a more holistic model of a cluster's MF. 

(i) Unresolved Binaries - It is well-established in the 
literature that many stars have equal or lower mass com- 
panions, with binary fractions ranging from ~ 0.35 to 1.0 
for massive O stars to ^0-3 to 0.5 for solar mass stars 
Preibisch et al.|1999l|Zinnecker fc Yorke|2007| |Sana| 



(e.g. 



et al.|201l[|2012[|Kiminki fc Kobulnicky|20l2] . SucITTF 
naries are typically difficult to detect, particularly with 
photometry, and are frequently treated as single stars. 
Several studies have shown that failing to account for 
unresolved binaries can cause systematic steepening to a 
cluster's MF slope by ~ 0.1 - 0.5 depending on the intrin- 



sic M F slope (e.g., Sagar fc Richtler|1991||Mafz Apellaniz| 
2008} 

Ln this paper, we have assumed that all data points are 
from single stars. However, there are at least two strate- 
gies for incorporating binaries into the presented frame- 
work. The first is to construct a model of the probability 
that a star is a binary based on measured luminosity and 
color, or estimated mass, of the primary star. It would 
be straightforward to employ such a model as a prior on 
the the mass PDF of each star. 

A second method would be to build binary star tracks 
into the process of photometric SED fitting, allowing one 
to compute the probability of a star being better de- 
scribed as a binary, in which case the mass PDF of the 
primary star would likely be changed dramatically, i.e., it 
would be flatter as one would expect little constraint on 
the mass of the primary solely from broadband photo- 
metric observations. The net effect of either approach 



is to increase the uncertainty in the mass of the pri- 
mary star, which reduces its contribution to constraining 
a cluster's MF. 

(it) Mass Segregation - Various observations have 
shown that massive stars are typically located in the cen- 
tral regions of young clusters, while lower mass stars are 
preferentially found in the outskirts of a cluster. The 
origins of observed mass segregation are often attributed 
cither to pr imordial star formation o r int ernal cluster dy- 



namics see 



Zinnecker fc Yorke|2007 and |Portegies Zwart 



et al. 2010 and references therein). The general effects 
ot mass segregation on MF measurements include flat- 
ter slopes in the inner regions, i.e., the stellar popula- 
tions have a higher percentage of high mass stars, and 
a steeper slope in the outskirts, where there are fewer 
high mass stars. Mass segregation can also introduce 
bias into the global MF measurement, depending on the 
treatment of completeness and the radius to which the 

data are considered. 

Inte restingly, simulations conducted by |Ascenso et al. 



(12009]) have shown that the effects of mass segregation on 
the MF slope are almost entirely degenerate with obser- 
vational completeness. That is, the completeness func- 
tion of a resolved cluster is a strong function of location 
such that the completeness limits are brighter in the cen- 
ter, and fainter in the outer regions. As a result, the 
recovery of a flatter MF slope in the central region of a 
cluster may simply be due to not properly correcting for 
a spatially variable completeness function such that the 
number of unobserved lower mass stars is significantly 
underestimated. 

With the presented probabilistic framework, the most 
straightforward way to account for mass segregation ef- 
fects in modeling a cluster's global MF is to ensure that 
the completeness function is well-characterized at all po- 
sitions within the cluster. Given that completeness mim- 
ics the behavior of mass segregation, an accurate ac- 
counting of spatial completeness will mitigate biases on 
the recovery of a cluster's global MF slope. At a funda- 
mental level, a well-characterized completeness function 
necessitates extensive artificial star tests, which can be 
computationally expensive. Within the PHAT program, 
large parts of the fundamen tal science require exte nsive 



such artificial star tests (see Dalcanton et al. |2012 ). As 



such, we anticipate having well-sampled spatially varying 
completeness functions, which will provide the necessary 
characterization of spatial variations in the cluster popu- 
lations to avoid strong biases in the recovered global MF 
slope. 

(Hi) Non-Coeval Populations - It is possible that young 
clusters are not purely single age, and instead have 
formed over some finite time interval. Applying a model 
that assumes coevality to a non-coeval population can 



slope (e. 




Miller & Scalo 


1979 


& Scalo 


2006), particularly whe 



The incorporation of an extended SFH into our proba- 
bilistic framework is fairly straightforward. It requires 
multiplying the intrinsic MF by an integral of a pa- 
ra meterized star formation history (cf. Equation 1 
in |Elmegreen fc Scalo] |2006[ ) , and then using MCMC 
sampling to simultaneously constrain the cluster's MF 
slope and star formation history, given the data. This 
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approach forms the basis of color-magnitude diagram 
model ing t echniques such as those di scussed in Dolphin 
( |2002| ) and |Harris fc Zaritsky| ( |2001| ), and is known to 
lead to very broad MF slope constr aints for obj ects with 
extended SFHs, e.g., galaxies (e.g., WeiszpOll ). 

(iv) Cluster Membership - In practice, there is always 
some degree of ambiguity as to whether a given star is a 
member of the cluster or whether it is a member of the 
surrounding field population. Thus, it is necessary to 
assign a membership probability to each star. Within 
the presented framework, the membership probability 
enters as a linear multiplication term that simply serves 
to weight each p(d | 9, obs), where the weighting factor 
ranges from zero (not a cluster member) to unity (def- 
initely a cluster member). The simplest way to assess 
the probability of cluster membership is through a sta- 
tistical comparison with a color-magnitude diagram of a 
surrounding field population. 

7. CONCLUSIONS 

We have presented a probabilistic approach to con- 
strain the MF of a young resolved stellar cluster. This 
framework allows the incorporation of uncertainties in 
the masses of individual stars that may arise from the 
conversion of observed fluxes into stellar masses, avoids 
binning the measurement in mass, explicitly deals with 
completeness functions that may not be a simple step 
function of stellar mass, and assigns meaningful error 
bars to the parameters of interest. 

Adopting a single-sloped power-law MF model, we ex- 
plored the ability of this approach to constrain MF pa- 
rameters using mock data. In particular, we found: 

• For highly idealized mock datasets (perfectly 
known masses and completeness), we recover the 
slope of the input IMF with no systematic biases, 
and to a precision that depends on the number 
of observed stars and the dynamic range of the 
observed masses, i.e., log(M maXi0bs /M comp ). We 
verified that one can only derive a lower limit 
estimate on M max , namely M maX! ob s , which is not 
particularly informative in terms of understanding 
the upper mass limit of the IMF. 



• We computed the theoretical precision, i.e., lower 
limit, to which a can be measured as a function 
of the observed mass range, log(M maXi0bs /Af comp ) 
and the number of stars, N, which we approxi- 
mate analytically as a 3rd order polynomial. We 
compared the theoretical precision with selected 
literature IMF studies and found that ~ 3/4 of 
the literature studies quoted error bars below the 
theoretical limit, usually by a factor of ~ 2. After 
assigning these studies new and larger error bars 
based on the theoretical limit, we computed the 
weighted mean and weighted l-a values from the 
literature studies and found (a) = 2.46 with a l-a 
dispersion of 0.35 dcx. The broad uncertainties 
indicate that the mean a is consistent with several 
common IMFs such as Salpeter, Kroupa, Kenni- 
cutt, and Scalo, suggesting that the current state 
of MF studies have little leverage on discerning 
between significantly different high mass IMF 



slopes. This finding reinforces the need for a 
large scale systematic study of the high mass IMF 
in order to provide the much needed empirical 
constraints. 



We then considered the case where the masses are 
not perfectly known. Specifically, we generated 
mock mass lists where each mass had a log-normal 
error distribution. Using the same datasets, we 
demonstrated the differences in MF recovery using 
models that did and did not account for the effects 
of uncertain masses. In general, we found that if 
all the mass uncertainties were small, the recovery 
of a was not substantially biased. However, in the 
case of intermediate to large mass uncertainties, 
a failure to model the mass uncertainties resulted 
in systematic biases in the slope of recovered MF, 
the magnitude and direction of which depend on 
the magnitude of the uncertainty, e.g., moderate 
uncertainties result in too steep of a slope, while 
large errors result in too flat of a slope. Similarly, 
not accounting for mass uncertainties led to severe 
systematic overestimates of M maX! i nput (or under- 
estimates of a cluster's true age). In contrast, 
applying models that account for uncertainties in 
the stellar masses resulted in ncar-pcrfcct recovery 
of the input IMF slope with precision that was 
comparable to the case of perfectly known masses, 
in all except the cases with the most extreme mass 
uncertainties. We found constraints on M max to 
be in agreement with those from the case of highly 
idealized data. 



• For a completeness function that grows linearly 
from to 1 within a certain mass range, we found 
that perfect knowledge of the completeness func- 
tion resulted in excellent recovery of the MF slope. 
However, both moderately and extremely incorrect 
completeness functions led to strong systematic 
biases in the recovered MF slope. We strongly 
recommend that extensive artificial stars tests be 
used in all resolved MF studies. In lieu of this 
option, only data that are nearly 100% complete 
should be utilized to minimize systematic biases, 
although the precision on the MF constraints will 
decrease in this case. 



• We discussed factors that can influence MF slope 
determination, but that were not included in this 
paper. Such effects include unresolved binaries, 
mass segregation, non-coevality, and cluster mem- 
bership. In each case, we suggest ways in which 
such uncertainties can be folded into the presented 
probabilistic framework. 
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APPENDIX 
A. MARKOV CHAIN MONTE CARLO 

Markov chain Monte Carlo (MCMC) algorithms are a broad class of numerical techniques for estimating the form 
of proba bility distributions. There are many textbooks ( |Bishop||2003| |Press et al. 2007 MacKay||2003| |Gelman"et 
al. 2003) that describe the general formalism and application of various MCMC methods. Here, we provide a brief 
overview of the general theory behind MCMC comment on the specific algorithm that we employ in this paper. We 
refer the interested reader to the above references for more detailed information. 

MCMC provides a prescription for drawing a set of unbiased samples from a probability distribution function (PDF) 
that can be evaluated (up to a normalization constant) given a set of parameters. In this paper, we use MCMC to 
estimate the distribution of model parameters (a, M, nax , and A pre( j) for a cluster's mass function (MF) that is consistent 
with the observations and marginalized over the nuisance parameters — M max and AT pre< j — when only the MF slope, 
a, is desired, for example. After running a MCMC chain and obtaining K samples 9^ — {ctk, M maXj fc, -/V pre d,fc}, the 
expectation value of a conditioned on the data and marginalized over 9_ a — {M max , A pre( j} is approximately given by 



a) = f ap(6\{Mi},N)d6- 



1 



N 
n=l 



(Al) 



Similarly, the marginalized PDF for a 



p(a\{M t },N)= p(9\{Mi},N)d9- a 



(A2) 

is given by the histogram of samples projected into the a plane. This marginalized PDF provides the desired constraint 
on the MF slope, while accounting for degeneracies with nuisance parameters. 

MCMC is generally much more efficient than alternatives such as rejection sampling or grid-based methods because 
it requires fewer calculations to provide a representative sampling of the density. 

The most common ly used MCMC algorithm is called the Metropolis-Hastings (M-H) algorithm (Metropolis et al. 
1953 Hastings 1970). While this is not the algorithm we use, it is instructive to outline how it works. The M-H 
algorithm "walks" around the parameter space stochastically starting from a user-defined initial position. Each step in 
the chain is determined by sampling a proposal position 9' from a (simple) distribution Q(9' | 0w) that only depends 
on the current position 9^ and then accepting (with replacement) this position with the probability 



A(9';9 {t) ) = mini 1, 



pie'iiMi^N) Q(eM\e>) 



p(eW\{Mi},N) Q(0'|0W) 



(A3) 



For this proje ct, we used an affine-invariant ensemble algorithm (emcea^ |Foreman-Mackey et al. 2012 Goodman 
& Weare 20101. This algorithm is more efficient than M-H when sampling a non-trivial density and it requires 
significantly less fine-tuning to achieve good performance. Instead of exploring the parameter space sampling a single 
point at a time, emcee includes an ensemble of L coupled "walkers" and the proposal distribution Q(9' e \ 9—g) for a 
particular walker I is based on the current positions of all the walkers. Specifically, an update step involves randomly 
choosing a walker from the complimentary set 9j 6 9_e = {9^ +1 \ • • • , 0\^\ 9f) 1 , • • • , 9^} and proposing a position 
along the vector between the two walkers 

9' = 9i + Z 



3 « 



where Z is a random variable distributed according to 



p(Z) <x 







if - < Z < a 
a 

otherwise 



for some a > 1. The proposal 9' t is then accepted (again with replacement) with the probability 



A{e' l ;ef ) ) = minjl. Z 
where D is the dimension of the parameter space. 



D _, p(9' e \{M l },N) 
p(9 { ; ) \{M l },N)^ 



(A4) 



(A5) 



(A6) 



ii 
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Table 1 

Literature Constraints on the IMF Slope 



Cluster log(AT) log(M max bs/M:omp) a a Min. Mass Max. Mass A °' i """' re Reference 





Name 










(M ) 


(M ) 






(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


1 


NGC2323 


3.31 


0.99 


2.94 


0.15 


0.4 


3.9 


1.66 


Kalirai et al. 2003 


2 


Mil 


3.26 


0.49 


2.49 


0.09 


1.1 


3.4 


0.32 


Santos ct al. 2005 


3 


NGC663 


3.15 


0.8 


2.38 


0.22 


1.6 


10.0 


1.54 


Pandcy ct al. 2005 


4 


NGC2168 


3.0 


0.78 


2.29 


0.27 


0.6 


3.6 


1.76 


Kalirai et al. 2003 


5 


NGC2422 


2.6 


0.44 


3.07 


0.08 


0.9 


2.5 


0.23 


Prisinzano ct al. 2003 


6 


Tr23 


2.48 


0.32 


2.4 


0.6 


1.0 


2.1 


1.42 


Bonatto et al. 2007 


7 


Stock2 


2.31 


0.44 


3.01 


0.4 


1.5 


4.1 


1.05 


Sanner et al. 2001 


8 


NGC4852 


2.3 


0.73 


2.3 


0.3 


0.6 


3.2 


1.38 


Carraro ct al. 2005 


9 


NGC654 


2.26 


0.98 


2.16 


0.05 


1.2 


11.5 


0.41 


Pandcy ct al. 2005 


10 


NGC4349 


2.21 


0.44 


2.18 


0.12 


1.8 


5.0 


0.3 


Tarrab et al. 1982 


11 


Orion 


2.13 


1.04 


2.39 


0.15 


1.3 


14.1 


1.31 


Tarrab et al. 1982 


12 


Pleiades 


2.1 


0.61 


2.99 


0.4 


1.0 


4.1 


1.31 


Sanner et al. 2001 


13 


Ly9 


2.09 


0.32 


2.1 


0.7 


1.1 


2.3 


1.41 


Bonatto et al. 2007 


14 


NGC3532 


2.06 


0.45 


2.55 


0.21 


1.6 


4.5 


0.5 


Tarrab ct al. 1982 


15 


Ly4 


2.02 


0.28 


2.3 


0.2 


1.0 


1.9 


0.37 


Bonatto et al. 2007 


16 


NGC5715 


1.97 


0.3 


2.3 


0.5 


1.1 


2.2 


0.92 


Bonatto et al. 2007 


17 


CygOB2 


1.97 


1.0 


1.9 


0.2 


10.0 


100.0 


1.48 


Massey et al. 1995 


18 


aPer 


1.95 


0.84 


2.73 


0.23 


1.3 


8.9 


1.14 


Tarrab et al. 1982 


19 


Pleiades 


1.94 


0.69 


2.89 


0.27 


1.3 


6.3 


0.93 


Tarrab et al. 1982 


20 


LH9 


1.92 


0.74 


2.4 


0.2 


10.0 


55.0 


0.78 


Massey et al. 1995 


21 


NGC346 


1.92 


0.85 


2.3 


0.1 


10.0 


70.0 


0.5 


Massey et al. 1995 


22 


Trl4/16 


1.91 


1.08 


2.0 


0.2 


10.0 


120.0 


1.72 


Massey et al. 1995 


23 


LH58 


1.82 


0.7 


2.4 


0.2 


10.0 


50.0 


0.66 


Massey et al. 1995 


24 


Praesepe 


1.81 


0.39 


3.12 


0.39 


1.3 


3.2 


0.72 


Tarrab et al. 1982 


25 


LH10 


1.81 


0.95 


2.1 


0.1 


10.0 


90.0 


0.6 


Massey et al. 1995 


26 


NGC2516 


1.81 


0.65 


2.95 


0.33 


2.0 


8.9 


0.97 


Tarrab et al. 1982 


27 


NGC6067 


1.8 


0.5 


3.18 


0.39 


4.0 


12.6 


0.85 


Tarrab et al. 1982 


28 


Haydes 


1.79 


0.33 


3.58 


0.48 


1.3 


2.8 


0.82 


Tarrab et al. 1982 


29 


NGC2264 


1.78 


0.56 


2.53 


0.29 


2.2 


7.9 


0.69 


Tarrab et al. 1982 


30 


NGC457 


1.71 


0.45 


1.36 


0.05 


5.0 


14.1 


0.09 


Tarrab et al. 1982 


31 


NGC6405 


1.68 


0.65 


2.84 


0.37 


2.0 


8.9 


0.99 


Tarrab et al. 1982 


32 


NGC2281 


1.67 


0.6 


4.17 


0.53 


1.4 


5.6 


1.27 


Tarrab ct al. 1982 


33 


NGC6281 


1.63 


0.45 


2.14 


0.23 


1.6 


4.5 


0.4 


Tarrab et al. 1982 


34 


NGC6633 


1.62 


0.29 


3.08 


0.45 


1.8 


3.5 


0.66 


Tarrab et al. 1982 


35 


Cz37 


1.61 


0.14 


-0.1 


0.9 


1.8 


2.5 


1.27 


Bonatto ct al. 2007 


36 


LH117/118 


1.6 


1.0 


2.6 


0.2 


10.0 


100.0 


1.08 


Massey et al. 1995 


37 


NGC2099 


1.59 


0.51 


2.6 


0.38 


2.2 


7.1 


0.71 


Tarrab et al. 1982 


38 


NGC4755 


1.58 


0.4 


2.1 


0.23 


5.6 


14.1 


0.36 


Tarrab et al. 1982 


39 


NGC6025 


1.57 


0.45 


3.06 


0.5 


2.5 


7.1 


0.84 


Tarrab et al. 1982 


40 


NGC2362 


1.56 


0.75 


2.3 


0.32 


2.5 


14.1 


0.96 


Tarrab et al. 1982 


41 


NGC2451 


1.54 


0.78 


2.1 


0.27 


1.3 


7.9 


0.86 


Tarrab ct al. 1982 


42 


IC2602 


1.53 


0.63 


1.88 


0.21 


1.3 


5.6 


0.48 


Tarrab et al. 1982 


43 


NGC7243 


1.48 


0.41 


1.83 


0.18 


2.2 


5.6 


0.26 


Tarrab et al. 1982 


44 


IC4665 


1.48 


0.54 


1.73 


0.16 


1.8 


6.3 


0.29 


Tarrab ct al. 1982 


45 


NGC6611 


1.48 


0.88 


1.9 


0.2 


10.0 


75.0 


0.71 


Massey et al. 1995 


46 


NGC1960 


1.4 


0.41 


1.24 


0.05 


3.5 


8.9 


0.07 


Tarrab ct al. 1982 


47 


NGC5662 


1.4 


0.6 


1.76 


0.19 


2.5 


10.0 


0.35 


Tarrab ct al. 1982 


48 


NGC2422 


1.38 


0.65 


2.38 


0.42 


2.0 


8.9 


0.84 


Tarrab et al. 1982 


49 


NGC1662 


1.38 


0.44 


2.16 


0.32 


1.8 


5.0 


0.44 


Tarrab ct al. 1982 


50 


NGC5460 


1.38 


0.35 


3.35 


0.71 


2.5 


5.6 


0.88 


Tarrab et al. 1982 


51 


IC1805 


1.38 


1.0 


2.3 


0.2 


10.0 


100.0 


0.72 


Massey et al. 1995 


52 


NGC2539 


1.36 


0.2 


2.61 


0.39 


2.0 


3.2 


0.44 


Tarrab et al. 1982 


53 


NGC7092 


1.36 


0.29 


1.72 


0.16 


1.8 


3.5 


0.19 


Tarrab et al. 1982 


54 


NGC3766 


1.36 


0.5 


1.53 


0.12 


4.5 


14.1 


0.18 


Tarrab et al. 1982 


55 


NGC2548 


1.36 


0.3 


4.05 


0.94 


2.0 


4.0 


1.11 


Tarrab et al. 1982 


56 


NGC6823 


1.36 


0.6 


2.3 


0.4 


10.0 


40.0 


0.71 


Massey et al. 1995 


57 


IC2391 


1.36 


0.65 


2.13 


0.34 


1.4 


6.3 


0.67 


Tarrab et al. 1982 


58 


NGC2301 


1.34 


0.6 


2.56 


0.49 


1.4 


5.6 


0.85 


Tarrab ct al. 1982 


59 


NGC884 


1.34 


0.5 


1.29 


0.06 


4.5 


14.1 


0.09 


Tarrab et al. 1982 


60 


NGC7160 


1.34 


0.8 


2.52 


0.45 


2.0 


12.6 


1.14 


Tarrab ct al. 1982 


61 


NGC4609 


1.32 


0.55 


2.02 


0.3 


2.8 


10.0 


0.46 


Tarrab ct al. 1982 


62 


Crl40 


1.32 


0.6 


2.58 


0.51 


2.0 


7.9 


0.85 


Tarrab et al. 1982 


63 


NGC2571 


1.3 


0.66 


3.3 


0.65 


2.2 


10.0 


1.19 


Tarrab et al. 1982 


64 


NGC2439 


1.3 


0.3 


1.6 


0.14 


7.1 


14.1 


0.15 


Tarrab et al. 1982 


65 


NGC1893 


1.28 


0.81 


2.6 


0.3 


10.0 


65.0 


0.7 


Massey et al. 1995 


66 


IC2581 


1.26 


0.64 


2.07 


0.36 


3.2 


14.1 


0.61 


Tarrab et al. 1982 


67 


NGC1528 


1.26 


0.45 


1.15 


0.03 


2.0 


5.6 


0.04 


Tarrab et al. 1982 


68 


NGC2482 


1.26 


0.51 


4.33 


0.95 


2.2 


7.1 


1.26 


Tarrab et al. 1982 


69 


NGC2251 


1.2 


0.49 


2.36 


0.49 


1.6 


5.0 


0.6 


Tarrab et al. 1982 


70 


NGC6242 


1.2 


0.5 


2.15 


0.39 


4.0 


12.6 


0.48 


Tarrab et al. 1982 


71 


NGC2232 


1.2 


0.6 


2.73 


0.63 


2.0 


7.9 


0.92 


Tarrab et al. 1982 


72 


NGC6531 


1.2 


0.76 


1.94 


0.34 


2.5 


14.5 


0.66 


Tarrab et al. 1982 


73 


NGC1664 


1.18 


0.26 


3.62 


0.96 


2.5 


4.5 


0.9 


Tarrab ct al. 1982 
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Table 1 — Continued 





Cluster 


log(JV) 


log(M max>obs /M comp ) 


a 




Min. Mass 


Max. Mass 




Reference 




(7 


Aa theory 




Name 










(M ) 


(M©) 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


74 


NGC581 


1.18 


0.35 


3.09 


0.78 


6.3 


14.1 


0.78 


Tarrab et al. 1982 


75 


NGC7790 


1.15 


0.45 


2.44 


0.54 


4.5 


12.6 


0.58 


Tarrab et al. 1982 


76 


NGC3590 


1.15 


0.4 


2.07 


0.36 


5.6 


14.1 


0.37 


Tarrab et al. 1982 


77 


NGC6709 


1.15 


0.2 


1.98 


0.28 


3,5 


5.6 


0.25 


Tarrab et al. 1982 


78 


NGC6871 


1.11 


0.35 


2.18 


0.4 


.5.6 


12.6 


0.38 


Tarrab et al. 1982 


79 


NGC5281 


1.11 


0.46 


2.02 


0.36 


3,5 


10.0 


0.38 


Tarrab et al. 1982 


80 


NGC2169 


1.08 


0.6 


1.74 


0.27 


2.8 


11.2 


0.35 


Tarrab et al. 1982 


81 


NGC1502 


1.08 


0,1 


1.56 


0.18 


5.6 


14.1 


0.17 


Tarrab et al. 1982 


82 


NGC1342 


1.08 


0.36 


3.46 


1.05 


2.2 


5.0 


0.96 


Tarrab et al. 1982 


83 


NGC2244 


1.08 


0.85 


1.8 


0.3 


10.0 


70.0 


0.55 


Massey et al. 1995 


81 


NGC6871 


1.04 


0.6 


1.9 


0,1 


10.0 


40.0 


0.51 


Massey et al. 1995 


85 


NGC7380 


1.04 


0.81 


2.7 


0.3 


10.0 


65.0 


0.51 


Massey et al. 1995 


86 


Berkeley86 


1.0 


0.6 


2.7 


0.4 


10.0 


40.0 


0.49 


Massey et al. 1995 


87 


NGC2323 


0.9 


0.2 


4.42 


1.73 


5.0 


7.9 


1.36 


Tarrab et al. 1982 


88 


NGC6913 


0.78 


0.6 


2.1 


0.6 


10.0 


40.0 


0.72 


Massey et al. 1995 


89 


CepOB5 


0.78 


0.48 


3.1 


0.6 


10.0 


30.0 


0.59 


Massey et al. 1995 



Note. — The literature IMF studies used in this paper. We have tabulated the number of stars (3) over the given mass 
range (4), and the reported values of a (5), where all values have been updated to reflect our usage of asaipeter = 2.35 and the 
1-ct uncertainty on a (6), allowing for a direct comparison with our calculated theoretical lower limits. The values in column 
(9) have been computed using the 1-a value listed in the literature, i.e., column (6) in this table, and the theoretical precision 
presented in Figure [8] There, we also show that nearly ~ 3/4 of the literature considered quote error bars smaller than the 
theoretical lower limit. 
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Table 2 

Polynomial Coefficients for A« Approximation 



i i 


Bi i 





2.90 


0.1 


-2.58 


0.2 


1.00 


3 


-0.25 


1.0 


-0.49 


1.1 


-3.51 


1.2 


5.00 


1,3 


-1.78 


2.0 


-0.05 


2.1 


1.30 


2.2 


-1.56 


2.; 3 


0.52 


3,0 


0.01 


3,1 


-0.09 


3.2 


0.09 


3,3 


-0.02 



Note. — Coefficients for the 3rd order polynomial approximation of t he theoretical precision on the IMF slope, Aa, as a function of 
the observed number of stars and observed mass range (see Equation |16[. 



